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ABSTRACT 

We apply a likelihood analysis to pulsar detections, pulsar upper limits, 
and diffuse background measurements from the OSSE and EGRET instruments 
on the Compton Gamma Ray Observatory to constrain the luminosity law for 
gamma-ray pulsars and some properties of the gamma-ray pulsar population. 
We find that the dependence of luminosity on spin period and dipole magnetic 
field is much steeper at OSSE than at EGRET energies (50—200 keV and 
>100 MeV, respectively), suggesting that different emission mechanisms are 
responsible for low- and high-energy gamma-ray emission. Incorporating a 
spin-down model and assuming a pulsar spatial distribution, we estimate the 
fraction of the Galactic gamma-ray background due to unidentified pulsars and 
find that pulsars may be an important component of the OSSE diffuse flux, but 
are most likely not important at EGRET energies. Using measurements of the 
diffuse background flux from these instruments, we are able to place constraints 
on the braking index, initial spin period, and magnetic field of the Galactic 
pulsar population, and are also able to constrain the pulsar birthrate to be 
between 1/(25 yr) and 1/(500 yr). Our results are based on a large gamma-ray 
beam, but they do not scale in a simple way with beam size. We estimate 
that about 20 of the 169 unidentified EGRET sources are probably gamma-ray 
pulsars. We use our model to predict the pulsar population that will be seen by 
future gamma-ray instruments and estimate that GLAST will detect roughly 
750 gamma-ray pulsars as steady sources, only 120 of which are currently known 
radio pulsars. 



1. Introduction 

Since the launch of the Compton Gamma Ray Observatory (CGRO) in 1991, pulsed 
gamma-rays have been detected from only 8 spin-driven pulsars. Because of the small 
number of detected pulsars, there remain many unanswered questions in pulsar gamma-ray 
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astronomy. The relationship between pulsar radio and gamma-ray emission is not clear. 
Aside from the Crab pulsar, the shapes of radio and gamma-ray profiles are quite 
different, with the peaks of the radio and gamma-ray profiles falling at different pulse 
phases, suggesting that low- and high-energy emissions come from different regions of the 
magnetosphere. In fact, Geminga, an extremely bright gamma-ray pulsar (GRP), is weak 
or quiet in the radio band. Furthermore, the gamma-ray profiles are broader than the radio 
profiles, suggesting that the beaming solid angle for gamma-rays is larger than that for 
radio waves. Models for the GRP emission mechanism, such as the polar cap and outer 
gap models, have been invoked to explain the features of GRPs, but no model has been 
successful at explaining all of the observations. 

Understanding GRPs may be the key to determining the nature of the unidentified 
EGRET sources, perhaps the largest mystery remaining from the CGRO mission. Of 
the 271 sources listed in the Third EGRET Catalog ( [Hartman et al. 1999 ), 169 are not 



identified as either pulsars or blazars, the only known EGRET discrete sources. While some 
of these sources are distributed isotropically and therefore expected to be AGN, 73 of the 
169 unidentified sources are located within 10° of the Galactic plane. Since pulsars are the 
only known Galactic EGRET sources, it is likely that these unidentified sources are pulsars 
or an as-yet-unrecognized source population. Modeling the GRP population is important 
for estimating how many of the unidentified sources we expect to be pulsars, which of the 
unidentified sources are most likely to be pulsars, and how we should go about detecting 
their pulsed emission. 

Understanding GRPs may also be important for explaining the intensity, Galactic 
distribution, and spectrum of the diffuse Galactic gamma-ray emission, the dominant 
feature of the gamma-ray sky. Accurate models for this diffuse emission are essential for 
detecting all but the very strongest gamma-ray sources. However, it is not clear whether 
this emission is due solely to cosmic ray interactions with interstellar material or partly 
due to unidentified point sources. Modeling the GRP population allows us to estimate the 
pulsar contribution to the diffuse emission. Furthermore, we can use measurements of the 
diffuse gamma-ray background to constrain important properties of GRPs. 

As the CGRO mission is ending and future gamma-ray missions are being planned, it 
is an appropriate time to assess the population characteristics of GRPs. Suitable modeling 
allows us to constrain luminosity laws, luminosity evolution, and spatial distributions of 
GRPs and to predict the pulsar population that future high-energy telescopes such as the 
Gamma Ray Large Area Space Telescope (GLAST), with an estimated launch date of 2005, 
will detect. 



In order to address the above issues within the confines of the small number of detected 
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GRPs, we have developed a likelihood analysis which employs all available information 
about GRPs, including pulsar detections, pulsar upper limits, and diffuse background 
measurements, to constrain their luminosity law and population parameters. In this paper, 
we describe the included data from the CGRO instruments OSSE (Oriented Scintillation 
Spectrometer Experiment) and EGRET (Energetic Gamma Ray Experiment Telescope). 
We present our model for pulsar gamma-ray luminosity and luminosity evolution, explain 
our likelihood analysis method, and discuss our results and their implications. 



2. The Data 
2.1. OSSE 

Three pulsars (the Crab and Vela pulsars, and PSR B1509— 58) have been detected 
by OSSE (|Ulmer et al. 1994| ; |Matz et al. ""15911 ; [Strickman et al. 199q ). Table 1 lists 



their 50—200 keV fluxes, along with upper and lower distance bounds, spin periods, period 
derivatives, surface dipole magnetic fields (i.e. B 12 = PP, for P in seconds and E>i 2 

in units of 10 12 Gauss), and total spin-down energy loss rate (i.e. E = IQQ, where / is 
the moment of inertia^ and Q is the spin frequency). We calculate the distance D to each 
pulsar using the Taylor and Cordes (1993) model for Galactic electron density, except for 
pulsars with more accurate distance measurements from parallax or a supernova remnant 
association. Following Taylor and Cordes, we include a 25% error in our distance estimation. 
We also include upper limits to the pulsed OSSE flux for 27 known pulsars ( [Schroeder et al. 



|1995| ; [Ulmer et al. 1999| ). We assume a beaming solid angle fi 7 of 2n, consistent with the 
uniformly wide pulse profiles observed for the detected pulsars. A photon spectrum with 
a power law index of —2 (i.e. F oc z/~ 2 ) has been assumed to convert all upper limits to 
the 50—200 keV energy range. While measured spectral indices for B1509-58 ( |Wilson et al. 



1993| ) and Vela ( |Strickman et al. 1996| ) are greater than —2, adopting a harder spectrum for 
the conversion of upper limits to luminosity changes the derived upper limits by only a few 
percent and does not affect the results of this analysis. We also make use of measurements 
of the OSSE diffuse gamma-ray background at 3 locations in the Galactic plane. Table 2 
lists the central longitude I of each of these pointings, the viewing periods used, and the 



measured diffuse flux ( Skibo et al. 1997 ; Kinzer et al. 1999 ). Viewing periods, or time 



intervals over which the spacecraft pointing was fixed, are defined in the Third EGRET 
Catalog ([Hartman et al. 1999|) . A power law spectrum with index —2 has been assumed to 
calculate the integrated 50—200 keV diffuse flux. Choosing a different spectral index does 
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not significantly change result sQ 



2.2. EGRET 



We include 7 EGRET pulsar detections in our analysis ([Nolan et al. 1996 ; Kaspi et~aT 



|1999|) . Table 3 lists the E > 100 MeV fluxes for these pulsars, along with distance bounds, 
periods, period derivatives, magnetic fields, and spin-down energies. We also include 353 
upper limits to the pulsed flux for EGRET pulsars (Nel et al. 1996; Kierro et al. 1995; 



Arzoumanian 199"7| ). Again, we assume a beaming solid angle of 2ir and a photon power law 



index of —2. Table 4 lists the included measurements of the EGRET diffuse gamma-ray 
background (|Hunter et al. 19971) for E > 100 MeV, calculated by assuming a spectrum with 



a power law index of —2. Choosing a different spectral index will change the derived diffuse 
fluxesQ the impact of this on the analysis is discussed further in Section 5.2. Individual 
viewing periods are not listed, as fluxes were calculated from an EGRET composite map. 



3. Luminosity and Population Model 



3 For a spectral index of —1.5, total diffuse flux over the 50— 200 keV band is 1.1 times that in Table 2. 
For a spectral index of —2.5, total diffuse flux is 0.9 times that given in Table 2. 

4 For a spectral index of —1.5, diffuse fluxes are 3 times that listed in Table 4. For a spectral index of 
—2.5, total diffuse fluxes are 0.5 times that listed in the table. 



Table 1. OSSE Pulsar Detections 



Name 


F 


A 


D u 


P 


P 




Bl2 


E 




(ergs kpc~ 2 s _1 ) 


(kpc) 


(kpc) 


(ms) 


(s/s) 




(10 12 G) 


(ergs/s) 


B0531+21 


3.33 ±0.4 x TO 34 


1.50 


2.50 


33.4 


4.21 x 10" 


-13 


3.8 


4.47 x 10 38 


B0833-45 


8.34 ± 5.6 x 10 31 


0.40 


0.60 


89.3 


1.25 x 10" 


-13 


3.4 


6.92 x 10 36 


B1509-58 


1.24 ±0.2 x TO 33 


3.50 


5.30 


150.6 


1.54 x 10" 


-12 


15.5 


1.78 x 10 37 
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Table 2. OSSE Diffuse Measurements at b = 0° 



1 


VPs 


F 


(degrees) 




(ergs kpc -2 s" 1 sr^ 1 ) 





5 + 16 


1.41 ± 0.2 x 10 37 


25 


7+13 


4.67+ 1.3 x 10 36 


95 


515 + 519 


3.53 + 1.2 x 10 36 



Table 3. EGRET Pulsar Detections 



Name 


F 


Di 


D u 


P 


P 




B\2 


E 




(ergs kpc~ 2 s _1 ) 


(kpc) 


(kpc) 


(ms) 


(s/s) 




(10 12 G) 


(ergs/s) 


B0531+21 


9.55 ± 0.96 x 10 33 


1.50 


2.50 


33.4 


4.21 x 10" 


-13 


3.80 


4.47 x 10 38 


J0633+1746 


3.53 ± 0.35 x 10 34 


0.12 


0.22 


237.1 


1.10 x 10" 


-14 


1.62 


3.24 x 10 34 


B0833-45 


6.78 ± 0.68 x 10 34 


0.40 


0.60 


89.3 


1.25 x 10- 


-13 


3.39 


6.92 x 10 36 


B1046-58 


2.39 + 0.57 x 10 33 


2.24 


3.73 


123.7 


9.59 x 10" 


-14 


3.47 


2.00 x 10 36 


B1055-52 


4.01 + 0.80 x 10 33 


1.15 


1.91 


197.1 


5.83 x 10- 


-15 


1.10 


3.02 x 10 34 


B 1706-44 


7.92 ± 1.59 x 10 33 


1.37 


2.28 


102.4 


9.30 x 10" 


-14 


3.09 


3.39 x 10 36 


B1951+32 


2.29 ± 0.46 x 10 33 


1.78 


2.96 


39.5 


5.84 x 10- 


-15 


0.49 


3.72 x 10 36 



Table 4. EGRET Diffuse Measurements at b = 



1 


F 


(degrees) 


(ergs kpc~ 2 s _1 sr _1 ) 


20 


1.70 + 0.30 x 10 41 


40 


1.44 + 0.25 x 10 41 


60 


7.85 ± 0.40 x 10 40 
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3.1. Luminosity and Period Model 

We model a pulsar's gamma-ray luminosity L as a power law in period and magnetic 

field, 

L = 1 p- a B 12 ", (1) 

where P is the period in seconds and B12 is the surface dipole magnetic field in units of 10 12 
Gauss. If the luminosity predicted by Eq. [I] is greater than some fraction e 7 of the pulsar's 
total spin-down energy E (i.e. E = IQQ oc PP~ 3 oc P~ 4 P 2 ), the model luminosity is set 
equal to e 7 P (see Section 5). 

Given a luminosity law and assuming a spindown law Q oc Q n , where Q is the spin 
frequency and n is the braking index, we may calculate a population-averaged gamma-ray 
luminosity, from which the total diffuse flux from GRPs may be derived. 

The spin period at age t of a pulsar born with initial spin period Po is 

P(t) = P (l+t/r ) 1/{n - l) , (2) 

where tq = Po/(n — l)Po is the initial spindown time. We assume the magnetic field remains 
constant throughout the pulsar's lifetime. Therefore, as B\ 2 oc \TPP for n = 3, P P = PP. 
However, while we expect n = 3 for magnetic dipole radiation, all measured pulsar braking 
indices are less than 3 QLyne et al. 1993"t [Kaspi et al. 1994| ; |Gouifles et al. 1992j ). If n does 



not equal 3, then the magnetic field estimate Bi 2 must be viewed as only an effective field 
strength that may differ substantially from the actual dipole field strength. 

For a given magnetic field strength B\2, the probability density function (PDF) for 
period P over a population of objects born at constant birth rate is 

f P (P\B 12 ) = (TgP)' 1 , P < P < P g , (3) 

where T g is the age of the Galaxy and P g = P (l + T 9 /r ) 1 ^ n_1 ' ) is the period of the oldest 
pulsar. We note that for fixed Bi 2 , P is uniquely determined by P. Throughout the 
analysis, we make the simplistic assumption that the PDF for Bi 2 is a delta function to 
keep the number of fitting parameters small. 

The total number of pulsars in the Galaxy is 

N psr = N psr T g = 10 8 iVo.oiT 9il o, (4) 

for a birth rate N psr = 0.01 yr^iVo.oi and Galactic age T g = 10 10 T 9 10 yr. The differential 
distribution of their spin periods is 

= N psr f P (P\B 12 ) = 2N psr T g PP- 2 = 10 8 - 3 s^AWT^oPP" 2 . (5) 
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Note that, because P g depends on n, this expression is dependent upon n. 



The luminosity PDF for a given magnetic field strength 5 12 , calculated from Eqs. [T] 
& i is 



f (T\n \ - -f p ( P \ Bl ^ _ 2 ^ P 2/9/a r -(l+2/o) 



and the population-averaged GRP luminosity is therefore 



(L\B 12 ) 



P 2 (a - 2) 



1 + 



T 



for a^2 and 1 



(6) 



(7) 



(9-2 p 2-a 



10 15 7 gl2 P -^ 

T,(a-2) 



for a > 2 and n = 3. 



3.2. Spatial Distribution 

We adopt a spatial distribution in order to calculate the total gamma-ray flux from 
pulsars in some direction through the Galaxy. We assume that pulsars are distributed in a 
Gaussian disk of radius r g with exponential scale height h. We also include a molecular ring 
of width w r at Galactic radius ry. We can express the pulsar number density at position x 
as a function of radius from the Galactic center r and height above the Galactic plane z as 

where rj is a dimensionless ring strength parameter and po = N psr /V g is the average Galactic 
pulsar number density, with N psr the number of pulsars in the Galaxy and V g the total 
Galactic volume. This expression is normalized according to / <i 3 xp(x) = N psr . We also 
define a dimensionless spatial distribution G(x) = p(x.)V g /N psr such that J d 3 xG(x) = V g . 
To convert the PDF in r and z (Eq. ||) to a PDF /d(-D) in distance from the Sun Z), we use 
the relation 

D = [r 2 + z 2 -2R Q rsm<p + R 2 Q ] l/2 (10) 

where r 2 = x 2 + for x = r cos and y = r sin 0, and i? Q is the distance of the Sun from 
the Galactic center. Assuming cylindrical symmetry, we integrate over 



fD(D) = 2^R- I S drdZp ^ z) \ 



>2 + z 2 + R 2 & _ D 2\ 1 

2R Q r 




where \r 2 + z 2 + R 2 & — D 2 \ < 2R Q r for all r and z. The integrals over r and £ are done 
numerically. 
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3.3. Flux Distribution 

The PDF for the flux F for a given gamma-ray luminosity L in terms of the distance 
PDF f D (D) is 

MF\L) = S^h = I ( 7^y /2 fn(D)F-^ (12) 



\dF/dD\ 2 \n 



7, 



for beaming solid angle fly. The flux PDF after integrating over a luminosity PDF, for 
fixed B12, is 

f F (F\B l2 ) = JdL f L (L\B 12 )f F (F\L). (13) 

Substituting Eq. ||, we do this integral numerically with integration limits determined 
by the luminosities corresponding to periods Pq and P g . If the luminosity predicted for 
period Pq is greater than e y E, the maximum luminosity is set equal to e 7 E. In this regime, 
fi(L\B 12 ) scales as L~ 3 ^ 2 instead of the L~( 1+2 / Q ) of Eq. |[ Including a maximum allowed 
gamma-ray efficiency is therefore important for all models with a < 4. 

The differential number of pulsars vs. flux, with beaming solid angle f2 7 , is 

= ^N psr f F (F\B 12 ), (14) 

which can be expressed as a histogram of F with logarithmic bins of width A log F 

dN, 



AN psr (F) = 2Fsinh(1.15AlogF)— (15) 

Eq. [15] gives the expected flux distribution for all Galactic pulsars. To calculate the 
expected flux from all pulsars in one beam (i.e. a cone of directions encompassing a small 
solid angle) centered on a particular direction, we take the average luminosity from Eq. [7] 
and integrate across the solid angle of the beam and our dimensionless spatial distribution. 
The total flux due to pulsars for a beam centered on Galactic longitude I and latitude b and 
of solid angle Qb is 

F(l, b, n b ) = J^(L\B 12 ) J dn b (l, b) J dD G(x), (16) 

where distance D from the Sun is related to r and z of Eq. || through 

x = -Dcos6sinZ (17) 
y = R Q — D cos 6 cos Z (18) 
z = Dsinb. (19) 
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4. Likelihood Analysis 

We use a likelihood analysis which exploits the detected GRPs, GRP upper limits, 
and GRP diffuse background measurements to find the best model to describe the GRP 
population. The total likelihood for model (i.e. one combination of parameters a, (3, 7, 
n, P , and Bu) is factorable into the likelihoods for detections, upper limits, and diffuse 
background measurements and can be expressed as 

£ tot (0) = £ dct (0)£ up (0)£ dif (0). (20) 

The total likelihoods for the detections, upper limits, and diffuse measurements are the 
products over the individual measurements such that 

N det 

£det(0) = II Atet,i(0), (21) 

8=1 

A IP (0) = II and (22) 

i=l 

N dif 

Aiif(0) = II 4s/,i(©)> ( 23 ) 

where Ndet, N up , and N^f are the number of included detections, upper limits, and diffuse 
pointings, respectively. 



4.1. Detected Pulsars 

To calculate (0), the likelihood of model given the ith pulsar detection, we use 
Eq. Q to calculate the predicted model luminosity L^Q) for the pulsar's measured P and 
P and for model parameters a, /3, and 7. For an assumed pulsar distance D, we calculate 
a predicted flux Fi(®,D) = L i (®)/fl~ f D 2 and compare this to the detected flux Fj. The 
form of the likelihood function is chosen to be a Gaussian with width determined by the 
flux measurement error api- Integrating over a distance PDF f(D) (n.b. different from the 
distance PDF of Eq. [Tl]) to account for uncertainties in distance estimation, we may express 
the PDF for model given detection i as 

CdetA®) = dD f(D) g((Fi(e, D) - a Fi ), (24) 

JDi 

where g(x,a) = (27rcr 2 ) _1 / 2 exp(— x 2 /2a 2 ). The likelihood for a particular model maximizes 
when the pulsar's predicted flux Fj(0, D) equals the detected flux Fj. We assume a uniform 
distance distribution (i.e. f(D) = 1/(D U — D{) for D t < D < D u and otherwise). 
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4.2. Upper Limits 

For each pulsar with a measured upper limit, we calculate the predicted luminosity 
given its measured P and P and for model parameters a, (3, and 7. For an assumed pulsar 
distance D, we calculate a model flux Fi(@,D) and compare this to the measured upper 
limit Fi. The form of the upper limit likelihood function is chosen to be flat with a one-sided 
Gaussian roll-off with width determined by the upper limit measurement error api- Again 
integrating over a distance PDF f(D), we express the PDF for model given upper limit 
measurement % as 

C up ,i{&) = / f(D) dD K ' — — 25 

Jdi J Fi dF g(F - Fi);a Fi ) 

for Fi(@,D) > Fi. For Fi(@,D) < Fi, the likelihood is constant, so that all models which 
predict fluxes less than the measured upper limit are equally likely. 



4.3. Diffuse Flux Measurements 

To determine the likelihood for model given the measured diffuse flux in the ith 
beam, we calculate the population-averaged pulsar luminosity for model parameters a, f3, 
7, n, Pq, and B 12 from Eq. [7], and then use Eq. [IB] to calculate -Fj(0), the total predicted 
flux from all pulsars in the beam. We compare this to the measured diffuse flux in the beam 
Fi, by allowing the pulsars to comprise a maximum fraction of the total measured diffuse 
flux. Choosing a flat function with a one-sided Gaussian roll-off with width determined by 
&Fi, the error on the measured diffuse flux in beam i, as our likelihood function, the PDF 
for model given diffuse measurement % is 

. g, e ,^g(f-^.);^,) 

~ f™F, dF S(F - fdF,);a Fi ) (26) 
for Fj(0) > e^Fj. For Fj(0) < e^i^, the likelihood is constant. 

The forms of the likelihood functions for detections, upper limits, and diffuse 
background measurements are illustrated in Figure 1. 



4.4. Grid Search and Marginalization 



To find the parameters which best describe the GRP population, we calculate the total 
likelihood (Eq. |20|) across a wide range of values for the parameters a, (3, 7, n, Pq, and B12 



and find that combination of parameters which maximizes the total likelihood. For each 
parameter, we calculate marginalized PDFs as 



Jexc 9i d® g(g) 

fd& £(©) 



(27) 



so that, for a given parameter 0^, the marginalized PDF is the normalized integral over all 
other parameters. A confidence interval C on parameter Qj can be calculated from Eq. |27] 
by finding the range of Qj in which a percentage C of the total probability is contained. 



We make several assumptions about the GRP population at the start of the analysis. 
For both the OSSE and the EGRET pulsars, we adopt a spatial distribution (as defined by 
Eq. |^) parametrized by r g = 6 kpc, r r = 4 kpc, w r = 1.5 kpc, and r\ = 0.25 (n.b. rj = 0.25 
corresponds to a 40% increase of pulsar number density at the radius of the molecular ring 
from its value for a Gaussian disk alone). For the OSSE pulsars, all of which have ages less 
than 10 4 years, we adopt a scale height of h — 0.1 kpc. For the EGRET pulsars, the oldest 
of which is 10 6 years old, we adopt a scale height of h = 0.5 kpc, as these older pulsars will 
have moved further away from their birth sites. We take R Q to be 8.5 kpc throughout. The 
corresponding PDF for distance D from the Sun, calculated from Eq. [11], is illustrated in 
Figure 2. We take the age of the Galaxy to be 10 10 years and take the pulsar birthrate to 
be 1 per 100 years. We set e 7 , the maximum fraction of spin-down luminosity a pulsar may 
convert into gamma-ray luminosity, to 0.5, consistent with the highest e 7 measured for the 
known pulsars (e 7 ~ 0.3 for PSR B1055-52). Finally, we set e<2, the maximum fraction of 
the diffuse flux attributable to pulsars, equal to 0.5. In Sections 5.1 and 5.2, we discuss the 
effects that changing these assumptions have on the analysis. 



Using the 3 OSSE detections, 26 upper limits, and 3 diffuse measurements described 
in Section 2.1, we calculate the total likelihood across a wide range of values for a, (3, 7, 
n, Pq, and Byi. The likelihood maximizes for one combination of parameters a, /3, and 7, 
and is a maximum across a range of values for n, P , and Bi 2 . Figure 3 shows contours of 
equal log likelihood as a function of a, (3, and log 7, taken in pairs. Figure 4 shows the 
marginalized PDFs for all 6 parameters. Best parameter values for a, /3, and log 7, with 
95% confidence intervals calculated from the marginalized PDFs, are 8.3 7.6t\i, and 



5. Analysis and Results 



5.1. 



OSSE Results 



-12- 



19.4l|g, respectively, leading to a best-fit luminosity law of 

L = 10 19A P- S - 3 B 12 7S ergs/s. (28) 

The marginalized PDFs allow us to place a 95% confidence upper limit on n of 4.75, a 
95% confidence lower limit on Pq of 6.5 ms, and a 95% confidence upper limit on B of 10 14 
Gauss. Changing e^, the fraction of the total diffuse flux that is allowed to be attributable 
to pulsars, from our assumed value of 0.5 alters these limits slightly. Allowing pulsars to 
comprise up to 100% of the total diffuse flux results in a lower limit of 4.5 ms on Pq, while 
allowing pulsars to comprise no more than 10% of the total diffuse flux leads to an upper 
limit of 3.75 on n. In Figure 5, we show how the fraction of OSSE diffuse flux due to pulsars 
varies with n, Pq, and Bi 2 . For best-fit values of a, /3, and 7, the population-averaged pulsar 
luminosity (Eq. [?]), and hence the diffuse flux from all pulsars, increases with increasing n, 
decreasing P , and increasing B 12 . Figure 5 shows that unidentified pulsars may be quite 
important in contributing to the OSSE diffuse background and should be considered in 
future models of diffuse Galactic emission. 

Given our best-fit luminosity law and assuming values for n, Po, and B\ 2 of 2.5, 10 ms, 
and 2.5 x 10 12 Gauss, respectively, we use Eq. [L^ to calculate the distribution of predicted 
OSSE flux for our model Galactic GRP population. According to this histogram, shown 
in Figure 6, and given OSSE's point source sensitivity of 6.3 x 10 31 ergs/s (in the 50—200 
keV energy range, and for an integration time of 5 x 10 5 seconds (|Gehrels, Chipman, & 



Kiiiffcii 1993|) ), OSSE could have detected 20 GRPs as steady sources in a full-sky survey. 



In Figure 6, we also plot the model-predicted OSSE fluxes for all known spin-driven pulsars 
with measured period and period derivative (includes 546 pulsars ( [laylor, Manchester, fcj 
Lyne 19913)). The pulsar with the fourth highest predicted flux in the OSSE energy range is 



PSR B0540— 69, the Crab-like pulsar in the Large Magellanic Cloud at a distance of 50 kpc. 
This illustrates the enormous capability a more sensitive low-energy gamma-ray instrument 
would have for detecting young, distant pulsars that are unlikely to be discovered through 
radio searches^. 

In Figure 7, we list the 50 spin-driven pulsars with highest predicted OSSE fluxes in 
order of decreasing predicted flux, along with measured fluxes or upper limits. Error bars 
are derived from the uncertainties on both the luminosity law parameters a, j3, and 7 and 
on the pulsar's distance. This plot allows us to predict which pulsars are good candidates 
to search for gamma-ray emission with more sensitive low-energy gamma-ray instruments. 
Because our analysis assumes that pulsars' beams point directly at us, a pulsar's actual 



5 We note that 10 hours of integration are necessary to detect B0540— 69 in the radio (Manchester et al. 
19931 ) . 
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detectable gamma-ray flux may be lower than that listed in Figure 7, depending on the 
orientation of the beam and our line of sight. However, if the successfully detected pulsars 
have beams that are somewhat skewed from our line of sight, then the actual flux of other 
pulsars may be larger than the prediction. 

5.2. EGRET Results 

Using the 7 EGRET pulsar detections, 353 upper limits, and 3 diffuse background 
measurements detailed in Section 2.2, we calculate the likelihood for a wide range of 
parameter values and again find a well-defined likelihood maximum in a, (3, and 7. Figure 
8 shows contours of log likelihood vs. pairs of these parameters. The best-fit luminosity law 
is 

L = 10 32 P- 1 - 8 B 12 15 ergs/s, (29) 

quite different from the OSSE best-fit law. Figure 9 shows the marginalized PDFs for all 6 
model parameters. Best parameter values for a, (3, and log 7, with 95% confidence intervals 
calculated from the marginalized PDFs, are 1.82^q'^, 1.54^q }5, and 32.041^, respectively. 
We can place an upper limit at 95% confidence on n of 3.6, a 95% confidence lower limit on 
Pq of 7.1 ms and a 95% confidence lower limit on £? 12 of 1.9 x 10 10 Gauss. 

In Figure 10, we show the fraction of EGRET diffuse flux due to pulsars for our 
best-fit model and for various values of n, P , and B\ 2 . Because the dependence of the 
population-averaged luminosity (Eq. [7]) on period and magnetic field is much different than 
for the OSSE pulsars, the total diffuse EGRET flux varies differently with these parameters 
and actually increases with increasing Pq and decreases with increasing B12 for n = 2.5. 
Figure 10 show that the contribution from unidentified pulsars to the diffuse flux is, at 
most, a few percent, and therefore is not important for models of the diffuse background. 
There are, however, several factors which influence this conclusion. For instance, as noted 
in Section 2.2, if the spectral index of the diffuse emission is —1.5 and not —2, integrated 
diffuse fluxes would be 3 times greater. In this case, pulsars might comprise as much as 
10% of the diffuse Galactic emission. Additionally, if the GRP scale height is lower than the 
0.5 kpc that we have assumed, a greater fraction of the diffuse flux could be due to pulsars. 
However, even with a scale height of 0.1 kpc, as used for the OSSE population, this fraction 
is always lower than 10% for reasonable values of n, Po, and B>i 2 . 

In Figure 11, we show the predicted flux distribution for a model population of GRPS 
given our best-fit luminosity law (Eq. [15]) and for values of n, Pq, and B i2 of 2.5, 15 ms, 
and 10 12 Gauss, respectively. Predicted fluxes for all known radio pulsars with measured P 
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and P are also shown. The point source sensitivity of EGRET (Se = 1 x 10~ 7 ph cm -2 
s~ 4 = 10 31 - 54 ergs kpc -2 s" 1 ) and the expected sensitivity of GLAST (S G = 4 x 1CT 9 ph 
cm -2 s -1 = lO 31,54 ergs kpc~ 2 s^ 1 ) are marked by dashed lines. We take Sq to be twice 
that calculated by Gehrels & Michelson (1999) for a 2-year all sky survey, recognizing that 
because pulsars are primarily a low-latitude population, surveys for them will suffer from the 
increased diffuse background in the Galactic plane. We estimate that EGRET should have 
detected roughly 20 pulsars. This number is tantalizing as, of the 73 unidentified EGRET 
sources within 10 degrees of the Galactic plane, 17 have non-variable fluxes, as we expect 
for pulsars (|McLaughlin et al. 1996| ; [McLaughlin fc Gordes 19991) . We therefore suggest, as 



have others, that these unidentified EGRET sources are GRPs. Unfortunately, due to the 
long integration times needed to detect only a small number of gamma-ray photons from 
these sources, searches for pulsars of unknown period in the EGRET data are prohibitive. 
Searching for radio pulsar counterparts will help to determine the nature of some of these 
sources, but because radio pulsar beams are evidently narrower than gamma-ray beams, 
many of these pulsars will be radio-quiet. Therefore, detecting these EGRET unidentified 
sources as pulsars will likely have to wait for a more sensitive high-energy instrument which 
will allow blind periodicity searches. 

According to Figure 11, GLAST should detect roughly 750 GRPs as steady sources. 
We estimate that 140 of these pulsars (including all unidentified EGRET sources, if they 
are indeed GRPs) could be detected through a blind periodicity search of GLAST dataQ 
We note that, using the same detection criterion, only 3 EGRET pulsars are detectable 
through a blind periodicity search of EGRET data. Figure 11 also shows that about 
120 of the approximately 1000 currently known radio pulsars should be seen by GLAST. 
With this large sample of radio and gamma-ray emitting pulsars, which will undoubtedly 
increase as the Parkes survey and other surveys near completion, one may compare radio 
and gamma-ray pulse profiles, enabling us to determine where in the magnetosphere the 
gamma-ray emission is coming from and providing important insights on the merits of 
outer gap and polar cap models. For the pulsars which are not detectable in the radio and 
are too faint to perform a blind search with GLAST data, x-ray observations may yield 
counterparts. As the angular resolution of GLAST will be better than that of EGRET, 
these searches should be much more feasible because fewer trial sky positions need to be 
considered in referencing photon arrival times to the solar system barycenter. 

Our prediction of 750 GLAST-detectable pulsars is dependent upon several 



6 We assume a one month integration, a duty cycle of 1/2, and require a minimum N^/Nt of 50, where 
N s is the number of pulsed source counts and Nt is the number of total counts, for detection, as in Mattox 
et al. (1996). 
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assumptions. We are taking e 7 , the maximum fraction of spin-down luminosity a pulsar 
may convert to gamma-ray luminosity, to be 0.5. If the actual value of e 7 is lower than 
this, fewer pulsars will be detectable by GLAST. For example, if e 7 is actually 0.1, we 
expect 610 GLAST-detectable pulsars. In Figure 11, we plot the predicted flux distribution 
for a population with e 7 of 0.01. While this value is so low as to be inconsistent with the 
measured gamma-ray luminosities of the known EGRET pulsars, the figure illustrates the 
effect a lower e 7 has on the expected flux distribution. We note that, because the scaling of 
the best-fit OSSE luminosity law with P and B\ 2 is far steeper than the scaling of E with 
these values, the predicted OSSE flux distribution is insensitive to e 7 . 

Our GLAST prediction is also dependent upon assuming a birthrate of 1/(100 yr). It 
is possible, however, that the actual pulsar birthrate is substantially higher or lower than 
this value. Adopting a much higher or lower birthrate while keeping all other parameters 
fixed would cause the histogram in Figure 11 to be inconsistent with the EGRET pulsar 
detections. However, because the flux histogram of Figure 11 depends not only on the 
number of pulsars but also on the average luminosity of a pulsar, changing the braking 
index, magnetic field, and initial spin period can compensate for a reduced birthrate. For 
reasonable values of these parameters (i.e. 2 < n < 3; 1 ms < Pq < 30 ms; 1 < B>i 2 < 5), 
birthrates as low as 1/ (500 yr) and as high as 1/ (25 yr) can be made consistent with the 
observed GRP flux distribution. 

Our histogram of predicted fluxes also depends on the spatial distribution model we 
choose for the pulsar population. Our model does not account for the location of the spiral 
arms or for local density variations, which is evidenced by the discrepancy between our 
model's predictions and the detections of the bright Vela, Geminga, and Crab pulsars. 
Adopting a model with a larger radial scale and/or including spiral arms would cause the 
number of GLAST-detectable pulsars and the expected density of pulsars in our region of 
the Galaxy to increase. 

Furthermore, and perhaps most significantly, because GLAST has not yet been built, 
its actual sensitivity is not yet known. Additionally, this sensitivity will be latitude 
dependent and will also depend upon models for the diffuse Galactic background. While 
we defer a more thorough analysis which includes these effects to a later paper, we do find 
that the number of detectable pulsars varies with GLAST sensitivity as Sq 1 . Therefore, if 
the actual sensitivity is half our assumed value, we expect roughly 1500 GLAST-detectable 
pulsars. For a sensitivity twice that we have assumed, we expect 500 GLAST-detectable 
pulsars. 

Figure 12 lists the known pulsars with the top 50 predicted EGRET fluxes, along 
with measured fluxes and upper limits, indicating which pulsars are good candidates for 
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gamma-ray detection with future instruments. Error bars are derived from the uncertainties 
on a, (3, 7, and distance. Again, as discussed in Section 5.1, a pulsar's actual detectable 
gamma-ray flux may be higher or lower than that listed in Figure 12, depending upon the 
correctness of our assumptions about beaming. We note that the measured and predicted 
fluxes for PSR B 1055-52 differ significantly. Some of this may be due to the uncertain 
alignment between the rotation axis and the line of sight to this pulsar ([Lyne fc Manchester 
1988| ; |Rankin 1993|) . There is also evidence that this pulsar is much closer than we have 



assumed. While we derive a distance of 1.5 kpc from its dispersion measure and the Taylor 
and Cordes (1993) model, Ogelman and Finley (1993), in their analysis of the ROSAT 
x-ray data, found that a distance of 500 pc would produce a more realistic estimate of the 
neutron star radius (15 km compared to 30 km for the 1.5 kpc estimate). Furthermore, 
Combi, Romero, and Azcarate (1997) derive a distance estimate of 700 pc from a study 
of the extended nonthermal radio source around the pulsar. Therefore, in Figure 12, we 
also mark the predicted flux for PSR B1055-52 for a revised distance of 500 pc. It is 
interesting that adopting a distance of 500 pc removes the constraint that this pulsar must 
be converting 30% of its spin-down energy into observable gamma-ray luminosity. With 
the revised distance estimate, a conversion percentage of only 3% is necessary. Similarly, 
we may also be overestimating the distance to Vela. While we are using the commonly 
quoted distance of 0.5 kpc for Vela ( [Taylor, Manchester, fc Lyne 1993|) , recent spectroscopic 
evidence suggests that this pulsar may actually be only 0.25 kpc away (Cha et al. 1999). 
Adopting this revised distance would alter our best-fit model modestly, but would have a 
large effect on the implied efficiency of gamma-ray production for the Vela pulsar, lowering 
it from 1.5% to 0.4%. Incidentally, assuming a 0.5 kpc distance for B1055-52, the Geminga 
pulsar is the most efficient converter of spin-down energy to gamma-rays, with an implied 
efficiency of 22%. The Crab pulsar, with an implied efficiency of 0.01%, is the least efficient 
of the EGRET pulsars. 



6. Discussion 

We have used OSSE and EGRET pulsar detections, measured upper limits, and 
diffuse background measurements to constrain the luminosity law and some properties of 
the gamma-ray pulsar population. The luminosity laws are quite different at OSSE and 
EGRET energies, suggesting that different mechanisms are responsible for pulsar low- and 
high-energy gamma-ray emission. The OSSE law [L oc P~ 8,3 Bi2 7S ) is roughly consistent 
with optical pulsar detections ( [Pacini fc Salvati 1983| ) and is steeper than that found for 
ROSAT/ASCA pulsar detections ( |Ogelman et al. 1993[ [Becker fc Trumper 1996[ |Saitc| 



|1997| ). It is quite different from models which predict that a pulsar's gamma-ray luminosity 
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should be some fraction of the total spin-down luminosity (i.e. L oc P^^B^ 2 ). The EGRET 
law (L oc P~ 18 Bi2 ' ) is more similar to this spin-down law and even more consistent with 
a luminosity law which scales with the open-field line voltage drop (i.e. L oc P~ 2 Bi 2 ). It is 
consistent with the basic polar cap ( Daughcrty fc Harding 1996| ) and outer gap ( Zhang £E 



Cheng 1997| ) models that have been published. More gamma-ray pulsar detections and a 
comparison between radio and gamma-ray pulsar beam shapes and properties are necessary 
to discriminate between models and the many versions of them. 

We note that a pulsar's age is proportional to P 2 Bu~ 2 , so that the gamma-ray 
efficiency, or the fraction of total spin-down energy converted to gamma-ray luminosity, 
decreases with increasing age for OSSE pulsars and increases with increasing age for 
EGRET pulsars. Because both the OSSE and EGRET best-fit luminosity laws differ from 
the standard E model, they produce a new ranking of pulsars as candidates for gamma-ray 
detection and result in new predictions for the pulsar population we expect to see with 
future more sensitive gamma-ray instruments. 

We find that a more sensitive low-energy gamma-ray instrument (50—200 keV energy 
range) would detect young, distant pulsars which are most likely not detectable through 
their radio emission. This would be enormously helpful for studies of young pulsars, pulsar 
populations, and supernova physics. A more sensitive high-energy gamma-ray instrument 
(0.1 — 100 GeV energy range) such as GLAST would detect roughly 750 GRPs as steady 
sources, most of which will have radio beams misaligned from our line of sight. Of these 
750 GRPs, roughly 120 will be known radio pulsars. About 140 of the 750 GRPs will be 
detectable as pulsars through blind periodicity searches, giving us a much more complete 
picture of the GRP population. GLAST will be able to do blind periodicity searches on 
all of the unidentified EGRET sources, of which 20 are probably GRPs. We note that all 
predictions are dependent upon an assumed e 7 and birthrate and will vary for different 
values of these parameters. While there is currently a large range of acceptable values for 
these parameters, comparing model predictions and actual GLAST detections will enable 
us to constrain these parameters much more tightly. Our prediction will also obviously 
depend upon the final sensitivity of GLAST, with the number of detectable pulsars scaling 
roughly as Sq 1 . 

We estimate the contribution of pulsars to the measured diffuse background in the 
OSSE and EGRET energy ranges, and find that while pulsars may be an important 
component of the diffuse background at OSSE energies, they are most likely not important 
at EGRET energies. We are able to use measurements of the diffuse background to 
constrain the braking index, initial spin period, and magnetic field of the pulsar population. 
The OSSE background allows us to place a 95% confidence upper limit of 4.75 on the 
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pulsar braking index, a lower limit of 6.5 ms on the pulsar initial spin period, and an upper 
limit of 10 14 Gauss on the pulsar magnetic field, while the EGRET background implies a 
95% confidence upper limit of 3.6 on n, a lower limit of 7.1 ms on P , and a lower limit of 
1.9 x 10 10 Gauss on B. We note that there are large errors on these parameters, and that 
the real distributions of these quantities are likely to be more complicated than we have 
assumed. Again, GLAST should allow us to constrain these parameters far more tightly 
than currently possible. 

In order to do a more sophisticated analysis of the GRP population, more GRP 
detections are necessary. Large strides in our understanding of GRPs will likely have to 
wait until the launch of GLAST. In the meantime, however, there are useful things that can 
be done to further our understanding. While searches for periodicities in EGRET data are 
very difficult, they may not be impossible if one takes an iterative approach (i.e. searching 
small, sequential intervals of data sets where fewer trial periods and positions are necessary) 
and uses Bayesian methods to detect pulses of unknown shape (e.g. [McLaughlin fc Cordes| 



1999| ). Searches for radio counterparts of the unidentified EGRET sources likely to be 
GRPs are also important, as some GRPs will emit in the radio and even non-detections 
provide useful information. Furthermore, whenever new x-ray or radio pulsars are detected, 
especially young or nearby pulsars, one should return to the OSSE and EGRET data and 
fold them. An accurate ephemeris may not be available but a known period allows a smaller 
range of parameters to search and a reduction of the search threshold. Some of the Parkes 
multibeam pulsars ( |Camilo et al. 19991) , for instance, are excellent gamma-ray candidates. 
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Fig. 1. — Likelihood functions for detections, upper limits, and diffuse measurements are 
plotted against predicted model flux, Fj. Dotted lines mark where predicted flux Fi equals 
the measured detection, upper limit, or diffuse measurement. 
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Fig. 2. — The PDF for distance D from the Sun for a model population of GRPs, calculated 
for a model parametrized by r g = 6 kpc, h — 0.1 kpc, r r = 4 kpc, w r — 1.5 kpc, and 77 = 0.25, 
and for R Q = 8.5 kpc. 
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Fig. 3. — Contours of equal log likelihood vs. pairs of parameters for the OSSE data. Crosses 
mark the maximum of the log likelihood. 
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Fig. 4. — Marginalized PDFs for all parameters for the OSSE data. 
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Fig. 5. — The fraction of OSSE diffuse flux attributable to unidentified pulsars as a function 
of braking index n, initial spin period Po, and surface magnetic field strength Bi 2 . One 
parameter is varied and other parameters are kept constant at nominal values of n = 2.5, 
Po — 15 ms, and B i2 = 1-0. The dashed line shows where the fraction of diffuse flux 
attributable to pulsars reaches unity and defines a range of parameter values which are 
disallowed. 
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Fig. 6. — Histogram of predicted OSSE pulsar flux for a = 8.3, (3 = 7.6, log 7 = 19.4, n = 
2.5, Po — 10 ms, and 5 12 = 2.5. The thick solid line shows the predicted flux distribution for 
a model population of GRPs. The dashed line denotes the point source sensitivity of OSSE. 
The thin solid line shows the predicted fluxes for known pulsars with measured P and P. 
Bins with highest predicted fluxes are labeled and correspond to the following pulsars: a) 
Crab; b) B1509-59; c) Vela; d) B0540-69; e) B1610-50, B1706-44. 
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Fig. 7. — Known radio pulsars listed in order of decreasing predicted OSSE flux. Crosses 
denote predicted gamma-ray fluxes, with error bars. Solid dots mark OSSE-detected fluxes, 
while open circles indicate upper limits. If the distance to B1055— 52 is indeed 500 pc, 
instead of the assumed 1.5 kpc, this pulsar's predicted flux would be 9 times higher, making 
it the 35th (as opposed to 47th) entry in this table. 
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Fig. 8. — Contours of equal log likelihood vs. pairs of parameters for EGRET data. Crosses 
mark the maximum of the log likelihood. 




Marginalized probability density functions for all parameters for EGRET data. 
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Fig. 10. — The fraction of EGRET diffuse flux attributable to unidentified pulsars as a 
function of braking index n, initial spin period Pq, and surface magnetic field B12. As one 
parameter is varied, other parameters are kept constant at values of n — 2.5, Po = 15 ms, 
and B 12 = 1.0. The dashed line shows where the fraction of diffuse flux attributable to 
pulsars reaches unity and defines a range of disallowed values for n. 
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Fig. 11. — Histogram of predicted EGRET pulsar flux for a = 1.8, (3 = 1.5, log 7 = 32.0, n 
= 2.5, Pq = 15 ms, and B12 = 1.0. The thickest solid line shows the predicted pulsar flux 
distribution of a model population of GRPs, assuming e 7 = 0.5. The thinner line illustrates 
the predicted distribution for e 7 = 0.01. Dashed lines show the sensitivities of EGRET 
and GLAST to steady sources and for a blind periodicity search. The thin solid line shows 
the predicted fluxes for known pulsars. Bins with highest predicted fluxes are labeled and 
correspond to the following pulsars: a) Vela; b) Geminga, Crab; c) J0631+105, B0656+14, 
B1509-58, B1706-44; d) B0950+08; e) B1046-58. 
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Fig. 12. — Known radio pulsars listed in order of decreasing predicted EGRET flux. Crosses 
denote predicted gamma-ray fluxes, with error bars. Solid dots mark detected fluxes, while 
open circles indicate upper limits. The triangle marks the predicted flux for PSR B1055-52 
for a distance of 500 pc. 



